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ABSTRACT 


A numerical investigation was conducted of viscous, compressible 
flow about a circular cylinder executing harmonic motion normal to the 
free stream flow direction. The oscillatory motion of the cylinder led 
to vortex shedding, and to a Karman vortex street. Calculations were 
made at a Mach number of .20, Reynolds numbers of 100 and 1000, a constant 
cylinder oscillation amplitude equal to 10% of the cylinder radius, and 
a range of cylinder oscillation frequencies. The variation with oscilla- 
tion frequency of the time-averaged-drag coefficient, root-mean- squared 
lift coefficient, and shedding frequency, are presented. In general, in- 
creasing the cylinder oscillation frequency decreased the time-averaged 
drag coefficient, increased the root-mean-squared lift coefficient, and 
had very little effect on the shedding frequency. 
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1.0 


INTRODUCTION 


1. 1 Background 

When launch vehicles are subjected to steady ground winds while 
in the vertical position, unsteady aerodynamic forces are exerted on them. 
These unsteady forces are associated with vortex shedding from the sides 
of the vehicle, and appear to be influenced by vehicle geometry, Reynolds 
Number, oscillatory vehicle motion, and surface conditions. 

In the past, specific launch vehicle configurations have been 
investigated by means of dynamically scaled wind tunnel models. However, 
wind tunnels are limited in size. As a result, no direct tests can be 
made of the aerodynamic characteristics of many vehicles of interest. 
Instead, scaled-down models of large vehicles are examined in wind tunnels, 
and the resulting data are extrapolated to full vehicle size on the basis 
of gasdynamic scaling laws. For example, by increasing the wind speed in 
a wind tunnel test of a scaled-down model, the Reynolds number is raised to 
levels characteristic of the full-sized vehicle. For a Newtonian in- 
compressible fluid, the two flows are then rigorously similar. However, 
wind speeds cannot be increased indefinitely without generating large 
perturbations to the flow field due to real-gas compressibility (in prac- 
tice such effects are to be expected for Mach numbers greater than .3). 
Since large scaling factors are needed for very large vehicles such as 
Saturn V, it has often been necessary in wind tunnel programs to violate 
the scaling laws. Typically, measurements are made at lower Reynolds 
numbers than desired, in order to avoid drastic changes to the flow field 
due to compressibility at high Mach numbers. A compromise is then struck 
in which wind tunnel tests are conducted at lower Reynolds numbers and 
higher Mach numbers than would characterize the true vehicle environment 
according to the gasdynamic scaling laws. In the case of the Saturn V, 
for example, wind tunnel tests have been conducted at Reynolds numbers of 

C. 

about 4.5 x 10 , while the corresponding full-scale Reynolds numbers should 
have been about 15 x 10^.'*' Practical difficulties of scaling wind tunnel 
test data for large vehicles provide a strong motivation for the develop- 
ment of alternate procedures for determining flow fields, such as numerical 
integration of the equations of gasdynamic motion. 
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To our knowledge, there are no adequate numerical or analytical 
methods for calculating the complete three-dimensional flow fields about 
such asymmetric vehicles as the Saturn V. However, it is realistic, using 
existing numerical tools, to attempt to calculate flow fields about full- 
scale launch vehicle configurations under conditions of plane two-dimensional 
flow. It is also feasible in such calculations to simulate the wind- 
induced motion of the vehicle by forcing the body cross section to oscillate. 
In earlier work, vortex formation and shedding were computed around station- 
ary body cross sections (right circular cylinders). The results were in 
quantitative agreement with experimental observation and with available 
independent theoretical calculations, to an accuracy which appeared to 
depend only on the length of machine time expended in computation. The 
numerical methods used for the computation were expected to be equally valid 
for flow around an oscillating body cross section. 
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1.2 


Objectives of the Program 

The primary objectives of this program were: 

a) To extend the numerical methods by which vortex flow was prev- 

2 

iously computed around stationary right circular cylinders, » to the 
case in which cylinders execute harmonic motion at right angles to the 
free stream flow direction. 

b) To calculate flow fields about oscillating cylinders. 

c) To determine the effects of cylinder oscillation and Reynolds 
number on the vortex shedding frequency and on aerodynamic forces around 
the cylinder. 
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1.3 


Technical Approach 


The calculations of vortex shedding about stationary cylinders, 
noted above, were effected using a computer code called "AFTON 2P" which 
integrates the general equations of continuum motion subject to the assump- 
tion of plane two-dimensional symmetry. The code is based on an "explicit" 
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finite difference representation of the equations of motion, ’ by 
which mechanical properties at a discrete set of points (known as a "finite 
difference mesh"), are updated from given initial conditions to any desired 
later stage of motion. For this program, AFTON 2P was again applied to 
the special case of time-dependent viscous compressible flow. 

1.3.1 The Problem of a Moving Interior Boundary 

In order to extend the AFTON 2P computer code to oscillating cyl- 
inder problems, the AFTON 2P equations had to be modified to include a 
closed Lagrangian boundary moving inside the region of flow. The motion 
of a closed Lagrangian surface in the interior of an Eulerian domain pre- 
sents a formidable numerical problem, to which two possible solutions were 
investigated. 

The first method considered was based on the unique capability of 
the AFTON 2P code to employ arbitrary time-dependent coordinate systems. 

Thus, it appeared feasible to define the surface of the cylinder as one 
coordinate surface of a time-dependent coordinate system, with the rest of 
the coordinate mesh topologically identical at any instant of time to the 
meshes used to calculate flows around stationary cylinders. The coordinate 
mesh near the cylinder would then have been almost stationary with respect 
to the cylinder, and would have included fixed points of the cylinder's 
surface. At the same time, coordinate lines emanating from these surface 
points would terminate at the space-fixed lateral boundaries of the system. 

To satisfy these conditions in a simple way the mesh points would all undergo 
periodic motion with the same period as the cylinder, but with amplitudes 
decreasing smoothly to zero at the boundary of the region of calculation. 
However, it was concluded after some study that the boundary value problem 
presented by the shuttling cylinder could be handled by a second method 
which was as rigorous as that just outlined, but which entailed much less 
computer programming effort. 
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The procedure adopted to accomodate cy linder motion took advantage of 
the fact that the system contained only one moving boundary; either the 
Eulerian flow region boundary could be considered fixed while the cylinder 
boundary moved, or vice versa. By making use of a coordinate system 
stationary with respect to the cylinder, the entire problem of coordinate 
motion was reduced (apart from the addition of a simple body force) to 
that of calculating transport of mass, etc., across the flow region boundary; 
variations from the conditions of free stream flow were relatively small at 
that boundary. In addition, by fixing the boundary of the region of cal- 
culation with respect to the cylinder, and not with respect to an external 
reference frame, logical problems due to the flow of mesh points in and out 
of the computational domain were avoided. Thus, all points of the finite 
difference mesh were made to undergo the same periodic displacement with 
respect to a fixed external reference frame as the center of the cylinder 
itself; all mesh points were tied to the cylinder in rigid-body fashion. 

While the AFTON 2P code had to be modified somewhat to provide it with a 
mesh of this kind, the code changes were minimal; the required modifications 
are discussed in the following sections. 

1.3.2 Body Fo rces in the Frame of the Oscillating Cylinder 

In the frame of the oscillating cylinder, the cylinder's acceler- 
ation is superposed on the flow field as an apparent or "kinematic" force, 
and a term appears in the momentum equation to account for this force. At 
any instant of time the resulting kinematic acceleration, which is constant 
over the whole flow field, is just the negative of the acceleration of the 
cylinder. 

The points of the cylinder were all displaced according to the 

equation: 


te = A Sin at (1) 

where Ax is the distance of displacement transverse to the direction of the 
free stream flow, and m and A are the frequency and amplitude of oscillation 
of the cylinder, respectively. The code was modified to calculate the 
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kinematic acceleration implied by Equation 1, and that acceleration was 
added to the acceleration terms previously taken into account in AFTON 2P. 

The momentum equation in the moving coordinate system is discussed in detail 
in Appendices A and B. 

1.3.3 Boundary Flow in the Frame of the Oscillating Cylinder 

The accelerated motion of the coordinate system also led to 
changes in boundary conditions at the perimeter of the region of calculation. 
Whereas frictionless sliding had been assumed at the lateral boundaries of 
the system in previous calculations of flow around a stationary cylinder, 
the fluid now had to be given a velocity component normal to those boundaries, 
equal and opposite to that of the cylinder in the external reference frame. 
Since this component of velocity varied periodically with time, finite dif- 
ference equations were formulated and programmed which describe material 
flow in either direction across each of the lateral boundaries. With 
fluid moving out of the region of calculation across a lateral boundary, 
that boundary became equivalent to the downstream boundary itself, and was 
described by downstream boundary equations developed earlier; with a minor 
modification for coordinate motion, exit conditions employed in this study 
were derived from the characteristic boundary condition used in previous 
stationary cylinder problems. Material entering the computational domain 
was simply assigned free stream properties (see Section 2.2 and Appendix C 
for a description of the lateral and downstream boundary condition routine 
employed in this program) . 

1.3.4 Conservation in the Frame of the Oscillating Cylinder; Summary of 

Code Changes Due to Coordinate System Motion 

In order to preserve the mass, momentum, and energy conservation 
properties of the AFTON 2P equations in the moving coordinate system, two 
additional modifications were required. First, transport of mass, downstream 
momentum, and energy across lateral boundaries were included in the finite 
difference equations for mass, momentum, and energy conservation; previously, 
AFTON 2P was programmed for these effects at upstream and/or downstream 
boundaries. Secondly, the work of kinematic body forces had to be included 
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along with the surface work of boundary stresses in an equation for total 
energy conservation. Details of the energy conservation calculation are 
presented in Appendices A and B. However, the cylinder velocity superposed 
throughout the region of flow had no effect on the calculation of kinetic 
energy, since the code in its original form placed no restriction on 
material velocity in the interior of the calculational domain. 

Thus, by viewing the flow from the frame of the cylinder, the 
equation and code changes necessitated by the cylinder's motion were a 
relatively slight generalization of the momentum and total energy calcula- 
tions, and the application at lateral boundaries of conditions already 
provided for upstream or downstream. 



1.4 


Summary of Results 


The principal results of the moving cylinder calculations lie in 
four areas. 

a) The asymmetries introduced by the oscillating cylinder were suf- 
ficient to cause vortex shedding without the introduction of artificial 
perturbations; such perturbations were required in previous stationary 
cylinder problems. ^ As part of this program, the effects of asymmetries 
due to transverse cylinder motion were studied at a Reynolds number of 
100, with particular reference to the shedding process. It was found that 
when the period of oscillation was much greater than the time, t Q , between 
the start of cylinder motion and shedding, shedding first occurred on the 
side of the cylinder which lay in the direction of initial transverse 
motion. For example, when the period of oscillation was equal to 3t OJ 
initial transverse motion of the cylinder in the positive x-direction (i.e., 
to the right with respect to the free stream) resulted in shedding of the 
right-hand vortex first. More generally, it was observed that regardless 

of the period of cylinder oscillation, vortex shedding first occurred on 
the side of the cylinder which accumulated the greater time-averaged shear 
stress in the interval t Q . 

b) On the basis of our numerical results for a Reynolds number of 
100, it was found that after periodic steady-state flow was achieved, the 
time-averaged drag coefficient, Cq, decreased slightly, and the Strouhal 
number, S, increased slightly, with increasing cylinder excitation frequency 
i (see Appendix D) . For instance, an increase in the dimensionless cylinder 
excitation frequency from 0.17 to 0.70, produced a drag coefficient decrease 
from 1.60 to 1.56, and a Strouhal number increase from S = 0. 17 to S = 0.18. 

c) Considerable root-mean-squared lift amplification took place when 
the dimensionless cylinder oscillation frequency, t>, was equal to or greater 
than the stationary cylinder Strouhal number, S. At a Reynolds number, R, 
of 100, the calculated value of S was .17, in agreement with experimental 
observation. With R = 100 and <k> = .17, a root-mean-squared lift coefficient 
of . 175 was calculated; the corresponding value computed for a stationary 
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cylinder « 0) is about .134. However, for = .70 and R = 100, the 
calculated root -mean-squared lift coefficient was 1.78 - more than ten 
times the corresponding stationary cylinder (j4 = 0) value. These results 
are discussed more fully in Section 2.6. 

d) On the basis of numerical results at = .17, it was found 

that the time-averaged drag coefficient, Cd, decreased with increasing 
Reynolds number. With fi = . 17 the value of Cp, was 1.60 for R = 100 and 
0.53 at R = 1000.0 - a result in qualitative agreement with experimental 
data for stationary cylinders. 


9 



2.0 


CALCULATIONS MADE AND RESULTS OBTAINED 


2. 1 Summary of Problems Run and. Description of Mesh Used 

Two stationary cylinder problems and three oscillating cylinder 
problems are considered in this report. The stationary cylinder problems 
were designated Problem 211.41 and Problem 131.0. The oscillating cylin- 
der problems were denoted Problem 133.0, Problem 133.1, and Problem 133.2. 

In Problem 211.41 the formation, development, and shedding of vortices from 
a stationary cylinder was calculated at R = 100.0. Some results of Problem 
211.41 were reported previously, ^ but lift and drag data are presented for 
the first time here. Problem 131.0 dealt with symmetric flow about a 
stationary cylinder and was run primarily to generate initial data for the 
moving cylinder problems. In Problems 133.0 and 133.2 the cylinder oscil- 
lation frequency was varied at R = 100.0. The Reynolds number was 1000.0 
in Problem 133.1, while the cylinder oscillation frequency was the same as 
that of Problem 133.0. In each problem the cylinder radius was 1.5 cm, the 
free stream Mach number was 0.20, and the amplitude of oscillation of the 
cylinder was equal to a tenth of its radius. A summary of parameter values 
for the problems of the program is presented in Table 1. 

The finite difference mesh which was used in all moving cylinder 
problems and in the stationary cylinder problem 131.0, consisted of 3520 
points and is shown in Figure 1. The upstream boundary is located 6.6 
diameters from the cylinder center, where it has little effect on the flow 
field in the neighborhood of the cylinder. The mesh was also used for 
Problem 211.41, except that the upstream boundary was located 2.25 diameters 
from the cylinder center. 
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2.2 


Boundary Conditions Employed 


In all calculations the density, specific internal energy, and 
velocity of material in the y-direction (the free stream direction; see 
Figure 1), were given their free stream values at the upstream boundary. 
For the stationary cylinder problems, the velocity of material in the x- 
direction was zero at the upstream boundary. Howeyer, since the moving 
cylinder problems of the program employed a mesh which was stationary in 
the frame of the cylinder (see Section 1.3), the x-component of material 
velocity at the upstream boundary was set equal and opposite to the cyl- 
inder velocity in the laboratory frame. 

As noted in Section 1.3.4, the method used to compute downstream 
boundary flow for all problems of this program was developed for previous 
stationary cylinder calculations, ’ and is based on the method of charac- 
teristics. The generalization to the case of interest here, in which the 
boundary moves with respect to the lab frame, is given in Appendix C. 

In the stationary cylinder problems, fluid was allowed to 
slide without friction along the lateral boundaries, i.e., at a lateral 
boundary the normal component of material velocity, and the tangential 
stress, were zero. For the moving cylinder problems, which were solved 
numerically in a frame stationary with respect to the cylinder (Section 
1.3), mass, etc., flowed across the lateral boundaries of the finite dif- 
ference mesh either from the free stream to the mesh, or from the mesh to 
the free stream. When the lateral boundary moved into the free stream, 
transport took place from the free stream to the mesh, and free stream 
conditions were imposed at the boundary. On the other hand, when the 
lateral boundary moved away from the free stream, then the characteristic 
boundary condition (Section 1.3.2) was used if mass flowed across the 
boundary out of the mesh, and free stream conditions were used if mass 
flowed into the mesh (see Appendix B for details of the characteristic 
boundary condition) . These rules led to stable numerical flow fields near 
the lateral boundaries of the oscillating cylinder problems (see Section 
2.5 for further discussion). 
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Finally, a no-slip boundary condition was imposed at the cylindrical 
surface in all problems. For the stationary cylinders, the velocity of fluid 
at the cylinder surface was then set equal to zero. For the moving cylinder 
problems, the surface velocity was zero in the computational frame, which was 
fixed with respect to the cylinder. 
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2.3 


Initial Conditions 


Apart from the no-slip condition imposed at the cylinder surface, the 
two stationary cylinder problems (Problems 211.41 and 131.0) were started from 
impulsive initial conditions, i.e., conditions of uniform velocity, density, and 
energy. The flow field about a full stationary cylinder at a Reynolds number 
of 100.0 (Problem 131.00; see Table 1), and at a characteristic time t of 2.7, 
was subsequently taken as the initial field of flow for the moving cylinder 
problems. At T = 2.7, a symmetrical vortex pair had developed in Problem 131.00, 
with a rear stagnation point located about 3.2 cylinder radii downstream of the 
cylinder center; the vortex centers were located about 1.855 cylinder radii from 
the cylinder center. The flow field for Problem 131.00 at a characteristic 
time of 2.7 is shown in the form of a velocity vector plot in Figure 2. The 
vectors of Figure 2 are proportional to the particle velocities at the points 
of the finite difference mesh; a mesh point is located at the tail of each 
vector. 
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2 . 4 Review of Numerical Results Obtained Previously for a 

Stationary Cylinder at a Reynolds Number of 100 

Vortex formation, development, and shedding were calculated about a 
stationary cylinder at a free stream Reynolds number of 100, and at a free 
stream Mach number of .20 (Problem 211.41; see Table 1). This system, which 
was otherwise symmetric about a plane through the cylinder axis, was perturbed 
asymmetrically to induce vortex shedding. The calculation was then carried 
through the shedding of five vortices. When a periodic state had been 
reached, the computed flow field contained a vortex sheet with a shedding 
frequency within 107, of the value measured for R = 100.0, namely, S = .17. 

The time-averaged drag coefficient, time-averaged lift coefficient, 
and root-mean-squared lift coefficient were also found from the stress field 
of Problem 211.41. For this purpose, the instantaneous force of the fluid 
on the cylinder was computed by integrating numerically the stress component 
normal to cylinder surface and the stress component tangential to the cylinder 
surface, over the cylinder surface. The component of the resultant force along 
the direction of free stream flow defines the drag force, and the component 
normal to the free stream flow direction is the lift force. The instantaneous 
lift and drag coefficients were then determined from the relations 


C 


D 


Drag Force 
%Po=U^ D 


( 2 ) 


Lift Force 

l D 


(3) 


where is the free stream density, is the free stream velocity, and D 
is the cylinder diameter. Finally, time-averaged lift and drag coefficients 
were found by numerical evaluation of the following integrals 
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dt 


( 4 ) 


C 


D 


°L = 



( 5 ) 


where the time t is measured from the start of cylinder motion. The root- 
mean-squared lift coefficient was determined by numerical evaluation of the 
integral 

C^Cr.m.s.) = 



For Problem 211.41, the time-averaged drag coefficient approached the 
value 0^ = 1.7, which is approximately 3% lower than the value of measured 
by C. Wieselsberger^. The time-averaged lift coefficient was zero; the instan- 
taneous lift coefficient oscillates about zero in a periodic manner. However, 
the calculated root-mean- square lift coefficient, C (r.m.s.), was .134. 
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2. 5 Vortex Street Calculation for an Oscillating Cylinder 

at a Reynolds Number of 100 

In the first moving cylinder problem attempted (Problem 133.0; see 
Table 1) the flow field was calculated about a right circular cylinder at a 
Reynolds number of 100.0. The cylinder was made to oscillate at a frequency 
corresponding to a Strouhal number of .17. As previously noted, the initial 
state of motion for Problem 133.0 was specified as the flow field of Problem 131.0 
at a characteristic time of 2.7; the cylinder displacement was prescribed by 
Equation (1) , and the double amplitude of oscillation was 107, of the cylinder 
diameter. Problem 133.0 was run to a characteristic time of 10.9, which corres- 
ponds to about 1.5 cycles of cylinder oscillation. As can be seen from 
Equation (1) , the cylinder moved initially in the positive x-direction (to the 
right), and shedding began with the right-hand vortex at a characteristic time 
of approximately 4.63 (about .85 msec after the start of the cylinder oscilla- 
tion). When shedding took place, (t = 4.63), the cylinder had completed about a 
third of its cycle of motion, with a time-averaged displacement. Ax, of .1085 cm 
in the positive x-direction. 

That vortex shedding took place first on the right side of the cylinder 
is consistent with the fact that the shear stress was larger on the right than 
on the left, in the period prior to shedding. A larger shear stress on the right 
than on the left is a physically reasonable outcome of the cylinder's initial left- 
to-right transverse motion; the right boundary layer is compressed while the left 
layer expands, and the normal gradient* of tangential velocity should therefore be 
larger on the right than on the left. It can be seen from Figure 3 that the 
numerical results bear out these observations; time histories of shear stress on the 
right- and left-hand cylindrical surface are presented in Figure 3, at points 
plus and minus ninety degrees from the forward stagnation point of the cylinder 
(i.e., at the points x = 1.5 cm, y = 0.0 cm, and x = 1.5 cm, y = 0.0 cm respec- 
tively). The time scale in Figure 3 is measured from the start of cylinder 
oscillation, when t = 1.19 msec and T = 2.7. From Figure 3 it is seen that 
the shear stress is greater on the right than on the left up to a time of .85 msec, 
or T « 4.1. 
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Vortex shedding in Problem 133.0 is illustrated by Figures 4-7, 
in the form of velocity vector plots. In Figure 4 (at a characteristic time 
of 7.08), the first right-hand vortex has moved far downstream, and the first 
left-hand vortex is also seen to be mbving downstream. The displacement of the 
cylinder at T = 7.08 is also exhibited in Figure 4, since the cylinder's center 
vras initially coincident with the indicated space-fixed origin of coordinates; 
velocity vectors and mesh points are all shown relative to the space-fixed frame. 
At a characteristic time of 7.95 (Figure 5), the first left-hand vortex is still 
propagating downstream, and a second right-hand vortex has begun to form. In 
Figure 6 (T = 8.9), the first left-hand vortex has moved still farther downstream, 
while the second right-hand vortex has become fully formed. The velocity vector 
plot of Figure 7 (T = 10.45) shows the second right-hand vortex shedding, and 
another left-hand vortex forming. As can be seen from Figures 4-7, the flow 
at the lateral boundaries is essentially free stream flow, with no evidence 
of numerical instabilities; the lateral boundary condition imposed in Problem 
133.0 therefore appears satisfactory for flows of this type. 

After many vortices have shed, the flow is experimentally observed to 
be almost periodic; a vortex detaches itself first from one side of the cylinder, 
then from the other side, and so on. Under conditions of near-periodic flow, 
times of arrival of successive vortex centers along a given line normal to the 
center line of the system, are almost equally spaced. 

In order to obtain a periodic flow field, Problem 133.0 was run far 
enough to permit the shedding of four vortices. From a careful examination of 
velocity vector plots, the trajectories of the centers of the four vortices 
were determined. These trajectories are shown in Figure 8, where y/a is plotted 
as a function of time("y" is the streamwise position coordinate and "a" the 
cylinder radius). The time between the passage of Left Hand Vortex One(LHVl) 
and Right Hand Vortex Two (RHV2) at y/a =3.9 was 1.26 msec, and the time period 
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between the passage of RHV2 and LHV2 at y/a = 1.6, was 1.28 msec (see Figure 8; 
shedding began on the right). Since these time periods are within 2% of each 
other, it is assumed that near-periodic flow exists 2.32 msec after the start 
of cylinder oscillation. By doubling the time differences cited, numerical 
approximations to the vortex shedding period were obtained; expressed as 
dimensionless Strouhal numbers (Appendix D) the period found at y/a =3.9 
was .175 and at y/a = 1.6 the period was .172. The measured Strouhal number 
for a stationary cylinder at R = 100.0 is .17. 

Based on Equations (2) through (6), values of the coefficients of 
time-averaged drag, time-averaged lift, and root-mean squared lift were computed 
for Problem 133.0. The time-averaged drag coefficient was 1.60 - about 6 % less 
than that calculated for a stationary cylinder at a Reynolds number of 100.0 
(see Section 2.4). The time-averaged lift coefficient was zero. The root-mean 
squared lift coefficient was .175, a value approximately 30% greater than the 
root-mean- squared lift coefficient calculated for a stationary cylinder at a 
Reynolds number of 100.0 (Problem 211.41; see Section 2.4). The lift amplifica- 
tion calculated for a cylinder oscillation frequency equal to the stationary 
cylinder shedding frequency, is consistent with available experimental results. 
Cincotta^ experimentally investigated unsteady aerodynamic lift forces on 
oscillating cylinders over a range of Reynolds numbers from .4 x 10^ to 10.7 x 
10 . Root-mean-squared lift amplication was observed when the model oscillation 
frequency was approximately equal to the stationary model shedding frequency. 
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2.6 


The Effect of Cylinder Oscillation Frequency on 
Vortex Shedding, Drag, and Lift 


To determine the effect of the cylinder oscillation frequency on the 
shedding frequency, time-averaged drag coefficient, and root-mean-squared 
lift coefficient, the frequency of oscillation of the cylinder was increased 
to an equivalent Strouhal number value of .70 (Problem 133.2, see Table 1). 

Initial conditions were again those of the flow field of Problem 131.0 at 
a characteristic time of 2.7 (see Section 2.3), and the Reynolds number was 
100.0. For Problem 133.2 vortex shedding is illustrated in the velocity 
vector plots of Figure 9 - 12; at the increased cylinder oscillation fre- 
quency of Problem 133.2, the velocities of the mesh points on the cylindrical 
surface are large enough to appear in the vector plots. 

Although the cylinder moved initially in the positive x-direction 
(to the right, as in Problem 133.0), the left-hand vortex ished first in Problem 
133.2, at a characteristic time T of 5.30. In the velocity vector plot of 
Figure 9 (T = 7.17), the left-hand vortex is displaced slightly downstream 
of the right-hand vortex. Shedding of the left-hand vortex first, despite the 
initial lef t-to-right cylinder motion, is consistent with the fact that the time- 
averaged shear stress was greater on the left than on the right at the time 
shedding took place. The possibility of such shear stress history arises when 
the transverse oscillation frequency exceeds the vortex shedding frequency. 

Other things being equal, if vortex growth is slight over a period of transverse 
oscillation, then circulation will accumulate on the left and right in about equal 
amounts in one cycle. Since the cylinder starts from its mean position, right to 
left motion begins after a quarter of a cycle is complete; circulation on the 
left will then tend to exceed that on the right during the third quarter-cycle, 
and any shedding in that interval would take place on the left. Actually, the 
initial incidence of shedding is complicated by vortex growth, and by other 
transients arising from both the initial conditions of flow and cylinder motion; 
circulation does not accumulate on the right during the first and fourth quarters 
of a cycle, in just the same way as on the left during the second and third quarters. 
As a result, to give a convincing account of initial shedding on the left may 
require a detailed understanding of relatively minor aspects of the motion, par- 
ticularly if the oscillation frequency is much greater than the shedding frequency. 
In Problem 133.2 the frequency of oscillation was about twice that for shedding. 
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Figure 13 shows shear stress histories for Problem 133.2 at points 
on the right- and left-hand cylindrical surfaces; the points are located ninety 
degrees (plus and minus) from the forward stagnation point of the cylinder (see 
Figure 1) . The time axis of Figure 13 is measured with respect to the start 
of cylinder oscillation. Based on the results presented in Figure 13, time- 
averaged shear stresses T (A) and T (r) were determined on the left and 
right sides of the cylinder, respectively , by numerical evaluation of the follow- 
ing integrals: 

-t 


T (A) 
xy 


t (r) 
xy 


o 

■if 


T (4) 

xy 7 


T (r) 
xy 


dt 


dt 


( 7 ) 


( 8 ) 


where t (A) and t (r) are instantaneous shear stresses at the left- and right- 
xy xy 

hand points, and t is the time. Shedding occurred 1.146 msec after the start 
of cylinder motion (i.e., at a characteristic time t equal to 5. 30), when the 
time-averaged shear stresses on the right and left, respectively, were 
T x y(r) = .00420 and T^(A) = .00455. In Problem 133.0 (see Figure 3) the 
instantaneous shear stress was always greater on the right prior to shedding, 
and the right hand vortex shed first (see Section 2.4). Thus, the results of 
Problems 133.0 and 133.2 are consistent with the idea that vortex shedding will 
occur on the side of the cylinder having the greater time-averaged shear stress 
at the time of shedding. 

Based on Equations (2) - (5) , time-averaged lift and drag coefficients 
were computed for Problem 133.2. The time-averaged lift coefficient was zero, 
i.e. , the instantaneous lift coefficient oscillates about zero. The time- 
averaged drag coefficient was 1.56 - about 11% less than that of a sta- 
tionary cylinder at a Reynolds number of 100.0. The time-averaged drag co- 
efficient is shown in Figure 14 as a function of cylinder oscillation frequency. 
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It can be seen from Figure 14 that the drag decreases with increasing oscillation 
frequency, a dependence that appears qualitatively correct. Due to the trans- 
verse velocity of the cylinder, which is directly proportional to the amplitude, 
A, and frequency of oscillation, w. (see Equation (1)), the instantaneous 
speed of free stream flow relative to the cylindrical surface will be equal to 
or greater than the free stream speed relative to a laboratory frame. Thus, 
the instantaneous Reynolds number, r' , relative to the cylinder, and based on 
the cylinder's diameter, will always be greater than or equal to the Reynolds 
number, R, relative to the laboratory frame. Therefore, if the time-averaged 
Reynolds number, R 1 , relative to the moving cylinder is defined as 



where 9 = U)t, then r' s r. Experimentally it has been found that the drag 

coefficient for a stationary cylinder varies inversely with Reynolds number*^. 
Therefore, the numerical result that at constant oscillation amplitude the 
drag coefficient decreases with increasing cylinder oscillation frequency 
appears at least qualitatively correct. 

The increase in the Reynolds number, R , with the cylinder oscillation 
frequency, u), can be used to establish a qualitative relation between u) and 
the time-averaged boundary layer thickness, 6, on the cylindrical surface. 

Since S' increases with the local Reynolds number per unit length around a 
cylindrical cross-section also increases as <a increases. Recalling that the 
boundary layer thickness, is inversely proportional to the square root of the 
local Reynolds number, it follows that 6 should decrease with increasing cylinder 
oscillation frequency. For an infinite flat plate executing simple harmonic 

g 

motion with frequency u) along the direction of free stream flow , the effective 
boundary layer at the plate has the thickness 
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where v is the kinematic viscosity; similar results apply for other solid 

surfaces oscillating in a fluid at rest or for a fixed surface in an oscilla- 

9 

ting stream . Thus, the decrease in boundary layer thickness with oscillation 
frequency, found in the results of Problems 133.0 and 133. 2, is not only 
physically understandable, but also has a firm theoretical basis in unsteady 
boundary layer theory. 

After a periodic steady-state was achieved, the vortices in Problem 
133.2 shed at a Strouhal number of .18. The variation of Strouhal number with 
cylinder oscillation frequency, for a constant double amplitude of oscillation 
equal to 10% of the cylinder diameter, is shown in Figure 15. As can be seen 
from Figure 15, increasing the cylinder oscillation frequency causes a slight 
increase in shedding frequency, as one might expect; the effective Reynolds 

_ i 

number R increases with frequency, and the Strouhal number is known experi- 
mentally to increase with Reynolds number. 

To establish a relation between transverse cylinder motion and the 
direction of the instantaneous lift force on the cylinder, the instantaneous 
lift coefficient was studied as a function of time for Problem 133.2 after a 
periodic state had been achieved. The time variation of the instantaneous 
lift coefficient, cylinder displacement, and cylinder acceleration are shown 
in Figure 16. A positive lift coefficient indicates a force in the positive 
x-direction (see Figure 1). Aside from a small phase angle, it can be seen 
from Figure 16 that the aerodynamic lift force on the cylinder generally 
opposes the cylinder acceleration, a result believed to hold for any cylinder 
oscillation frequency. 

From the instantaneous coefficients of Figure 16 and Equations (3) and 
(6), a root-mean-squared lift coefficient of 1.78 was computed for a cylinder 
oscillation frequency corresponding to a Strouhal .number of 0.70 (Problem 133.2). 
The root-mean- squared lift is therefore about ten times larger for a cylinder 
oscillation frequency corresponding to a Strouhal number of 0.70 (Problem 133.2), 
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than for an equivalent Strouhal number of 0.17 (Problem 133.0) --a tremendous 
lift increase. The rapid growth of lift with oscillation frequency reflects 
a corresponding growth in the transverse pressure difference across the cylinder, 
as the oscillation frequency increases. The pressure distribution on the 
cylindrical surface, at a cylinder oscillation frequency corresponding to 
>5= 0.17 (Problem 133.0), and at a time when a periodic motion was achieved 
(T = 10.45; see Figure 7), is presented in Figure 17. In Figure 17 the pressure 
is greater on the right-hand cylindrical surface than on the left-hand surface, and 
the lift force is therefore in the positive x-direction. The area between the 
curves defines the pressure force contribution to the lift and corresponds 
to a positive pressure lift coefficient, C equal to 0.127, a value approxi- 

J— i XT 

mately 24 percent less than the total instantaneous lift. The pressure distri- 
bution on the cylindrical surface at a cylinder oscillation frequency corres- 
ponding to 6 = 0.70 (Problem 133.2), and at a characteristic time T of 10.05, 
is shown in Figure 18. A velocity vector plot of the flow field at this tirie 
is presented in Figure 12. In this case the pressure is generally higher on 
the right-hand surface than on the left, and the lift force is in the negative 
x-direction. The area between the right- and left-hand curves now corresponds 
to Cj^p = - 0.804, a value which is about 23 percent smaller than the total 
instantaneous lift, but about 6.3 times that cited above for Problem 133.0. Such 
a result is qualitatively reasonable; as the cylinder oscillation frequency in- 
creases, the maximum compression of fluid increases on the side of the cylinder 
in the direction of motion, while the minimum compression decreases on the 
other side. Lift forces should therefore grow with cylinder oscillation fre- 
quency, as the numerical results indicate. 

The variation of the root-mean- squared lift coef f icient ,C^ (r .m. s . ) , with 
cylinder oscillation frequency is presented in Figure 19. The curve faired 
through the data is based on the experimental results of Cincotta'*', which show a 
peak value of C (r.ra.s.) at the aerodynamic shedding frequency for a stationary 
cylinder (S = 0.17); C^(r.m.s.) then decreases with S reaching the stationary 
cylinder root-mean-squared lift coefficient at approximately twice the aerodynamic 
Strouhal number for a stationary cylinder. Cincotta's data are also shown in the 
figure. 
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2.7 


Effects of Reynolds Number Variation For an Oscillating Cylinder 


To investigate the influence of Reynolds number on the shedding 
frequency, drag, and root-mean-squared lift, the flow field about an oscilla- 
ting cylinder was calculated at a Reynolds number of 1000.0 (Problem 133.1). 

The cylinder oscillation frequency and amplitude were the same as those of 
Problem 133.0. The finite difference mesh consisted of the two-dimensional 
array of 44 x 80 points shown in Figure 1. The initial flow field was again 
identical to that of Problem 131.0 at a characteristic time of 2.7, when a 
symmetrical vortex pair had already developed. The sudden increase in 
Reynolds number caused apparently random oscillations in the pressure on the 
cylindrical surface, and subsequently in the instantaneous drag coefficients. 

After approximately 250 timesteps of calculation these oscillations appeared 
to have damped out, and the flow was assumed to be that appropriate to a 
Reynolds number of 1000.0. The calculation was continued for 1.5 shedding 
periods. 

Vortex shedding in Problem 133.1 is illustrated by the velocity vector 
plots of Figures 20 - 23. As in the case of Problem 133.0, the right-hand 
vortex shed first, but now at a characteristic time t of 4.63. Initial shedding 
on the right was due to initial transverse motion of the cylinder in the 
positive x-direction with a period much longer than the time needed to induce 
shedding (see Section 2.6). A comparison of the wake at a Reynolds number of 
1000.0 (Figure 23, T = 9.08) with that at a Reynolds number of 100.0 (Figure 7, 

T = 10.45), indicates a more random set of velocity vectors for R = 1000.0. 

The turbulent character of the wake at R = 1000.0 is consistent with our previous 

3 

results regarding pressure oscillations recorded on the leeward side of a cylinder , 
and with the experimental results of Roshko"*. Earlier numerical calculations 
resulted in pressure variations which were smooth at R = 100.0, but exhibited 
elements of randomness at R = 1000.0; Roshko found that the vortex street from 
a cylinder was laminar at a Reynolds number of 100, and turbulent at a Reynolds 
number of 1000.0. 
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Based on the numerical results of Problem 133.1, the effects of 
Reynolds number on the shedding frequency, drag, and lift were determined. 

When a periodic steady state had been achieved, the vortices shed at a 
Strouhal number of 0.21, which corresponds to the Strouhal number measured 
for a stationary cylinder at R = 1000.0. The time-averaged drag coefficient 
was 0.53. A drag coefficient of 0.53 is approximately half that of a sta- 
tionary cylinder at R = 1000.0; this result is consistent with the drag 
dependence found at R = 100 as the cylinder oscillation frequency increased. 

The time-averaged drag coefficient is presented in Figure 24 as a function of 
Reynolds number for a cylinder oscillation frequency corresponding to a 
Strouhal number of 0.17. The stationary cylinder data presented in Figure 24 
were used as a basis for fairing the curve shown for the oscillating cylinder. 
It can be seen from Figure 24 that, as in the case of a stationary cylinder, 
the drag coefficient decreases with Reynolds number. In the oscillatory case, 
the computed root-mean-squared lift coefficient was about 0.18; the calculation 
for a stationary cylinder at R = 1000.0 was part of an earlier program in which 

C (r.m.s.) was not computed, nor to our knowledge have root-mean- squared lift 
L 

data been measured. 
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3.0 CONCLUSIONS AND RECOMMENDATIONS 


The work summarized in this report constitutes our first attempt 
to extend the AFTON 2P code, and the numerical methods it embodies, 
to problems of viscous, compressible flow about accelerating structures 
The principal conclusion reached in the work reported here is similar 
to that of our stationary cylinder calculations^; up to a Reynolds 
number of 5000 complete flow fields about accelerating bodies can, at 
feasible cost, be predicted quantitatively by AFTON 2P with an accuracy 
sufficient for almost any practical purpose. Furthermore, by incorpor- 
ating a variable timestep into the numerical technique, it is quite 
likely that an accuracy comparable to that obtained in the problems 
reported here can be achieved up to a Reynolds number of 15 x 10^ 

(the Reynolds number for a full scale Saturn V vehicle) and for a 
class of two-dimensional accelerating bodies which includes most shapes 
of practical interest. 

These conclusions rest both on numerical results obtained previous 

2 3 

ly for stationary cylinders », and on the specific results discussed 
in detail in Section 2.0. The basic supporting facts will now be 
summarized. 

Our previous stationary cylinder calculations were performed in 
two separate programs. In the first program, it was established that 
for right cylinders of fairly simple shape, important gross para- 
meters such as lift, drag, and Strouhal number could be calculated up 
to Reynolds numbers of 1000.0. For example, at a Reynolds number of 
100, the calculated drag coefficient was within 3% of the experimental 
value and Strouhal number was within 10% of the measured value (see 
section 2.4) . However, no significant effort was devoted to determin- 
ing the accuracy of the numerical predictions for the pressure and 
velocity fields. In the second program, detailed calculations were 
performed at Reynolds numbers of 100 and 1000, with finite difference 
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meshes that were much finer than any used in the first program. The AFTON 
boundary layer profiles along normals to the cylinder were nowhere in 
error by as much as 1%, as determined by comparison with the known exact 
solution. The pressure coefficient distribution around the cylinder at 
a stage of substantial vortex cevelopment, for a Reynolds number of 100 
and a Mach number of .2, was within 15% of an independent incompressible 
viscous flow calculation; the AFTON 2P fluid is compressible. Furthermore, 
a study of the numerical solution error showed that the mesh point density 
could be made sufficiently large in practice to reduce the maximum pressure 
coefficient error for the entire cylindrical surface to no more than 1.4%. 
Thus by employment of finer meshes, any desired solution accuracy can be 
obtained. 

In the oscillating cylinder problems presented in Section 2.0, 
vortex shedding was initiated by the asymmetries introduced by the shut- 
tling cylinder, without the introduction of artificial perturbations. 

This result conforms to experimental observation. The finite difference 
mesh employed in these oscillating cylinder problems was identical to 
that employed in the stationary cylinder problems of our first program. 
Thus, it is reasonable to assume that the accuracy of the calculated lift, 
drag, and Strouhal number in the oscillating cylinder problems is appro- 
priate to the corresponding stationary cylinder parameters calculated in 
the first program. Furthermore, the calculated relationships between the 
cylinder oscillation frequency and the root-mean-squared lift coefficient, 
time-averaged drag coefficient and Strouhal number are qualitatively con- 
sistent with clear and simple physical considerations (see Section 2.6). 
Finally, since increased accuracy was achieved at feasible cost by mesh 
refinement in the stationary cylinder problems of our second program, it 
appears that a similar increase in accuracy can be attained by employing 
finer meshes in the moving cylinder problems. 
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A development important to future calculations lies in the discovery 
of finite difference equations which preserve all the properties of 
self-consistency of the present equations, but in which the time-step 
can be varied from zone to zone. In this way, the variables of motion 
need only be calculated for a given zone as often as required by con- 
siderations of stability for that zone and its immediate neighbors, 
rather than for the least stable zone of the finite difference mesh. 

The increase in speed of solution of problems like those reported here 
was potentially so great that variable timestep (VTS) equations were 
investigated in one space dimension. An increased speed of solution 
by a factor of five was obtained. This factor should be even greater 
in two dimensions, using a two-dimensional counterpart of our one- 
dimensional VTS scheme. It is therefore recommended that a variable 
timestep be incorporated into the AFTON 2P computer code for moving 
cylinders, and that a numerical investigation be conducted of flow 
fields about oscillating cylinders in the higher Reynolds number regime 
(i.e., 10^ to 15 x 10 6 ) . The experimental data of J.J. Cincotta could 
be used to establish the accuracy of the numerical results in this 
range of Reynolds number. 
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APPENDIX A 


EQUATIONS OF MOTION IN A MOVING COORDINATE SYSTEM 

In this section the integral equations of motion are developed 
in a coordinate system fixed relative to a moving obstacle immersed 
in a viscous, compressible fluid. 

The integral equations of motion in the moving coordinate 
system will be derived in Cartesian tensor notation^, where the 
position vector is defined by x^ and the velocity vector is defined 
by u. (i = 1, 2, or 3). Consider a volume t of three-dimensional 
space having a surface area P which is a fixed relative to a labor- 
atory frame. In partial differential form, the conservation equations 
are as follows relative to the laboratory frame: 

Continuity 


t! + sir (p v 


0 


( 11 ) 


where t is the time, p is the density, and, by the summation conven- 
tion, the repeated indices indicate the summation 


+ r=- (pu,) 
3 


3^ < pu *> = 3^ (pu l> + 3^ < PU 2 ) ' ax, VK “3 

First Law 


3 (pe) , 5(pu P^ = p 

& t 9xp Pa dxg 


( 12 ) 


where E is the specific internal energy and P^ is the stress tensor. 
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Momentum 


d(pu ) 
v a' 

at 


a 

Sx Q 


( pu pV 



(13) 


The transformation equations relating the laboratory and 
moving frames will now be presented. Consider a moving 
coordinate system with a position vector x^ and velocity 
vector u. measured relative to this system. Let the velocity of 
the moving coordinate system be denoted V.(t) (i = 1, 2, or 3). The 
relationship between the laboratory frame and the moving frame can 
be defined by the following Galilean transformation: 



(14) 


where t' is the time in the moving frame. Based on the transformation 
equations ( 14 ) the partial derivatives in the laboratory and moving 
frames may be related as follows: 



3 

dt' 


a 

at- 



(15) 


Finally the particle velocities in the laboratory and moving frames 
are related as follows: 


u. = u. + V. (t) 
i l l 


(16) 
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Based on the above transformation equations, the conservation laws 
relative to the laboratory frame Equations (JL1 ) , (-12'.) and (13), 
can be transformed to the moving frame. The conservation equations 
relative to the moving frame are as follows: 

Continuity 

II. + sij «S> - o <«> 

First Law 


d (PE) _a_ 

at > a X £ 1 


(pEug) 


afu' 

a 


r to 3x 


3 


(18) 


where the 

~S^ 


stress tensor remains unchanged since it is proportional 
which is unaffected by the transformation. 


Momentum 


9P - U a 


TF + 4' (Pu 8 u a ) + p dF~ - tofT 

It is seen from the transformed Equations (17), (18) and (19) that, 
as discussed in Section 1.3, only the momentum equation has changed 
in the moving coordinate system. From the Gauss Theorem, relating a 
volume integral to a surface integral, Equations (17), (18) and (19) 
can be transformed to integral form. The Gauss Theorem is as follows: 


5l. 

dt 


3P, 


to 


(19) 



where T' is the fixed volume T defined relative to the moving coordinate 
surface (T 1 = T) , F' is the surface area of F defined relative to the 
moving coordinate system (r'=F), g is any tensor defined within T' 
and whose derivatives are continuous there, and is the direction 
cosine of the normal to the surface T ' in the direction P. Based on 
Equation (20) the final integral equations of motion relative to 
moving coordinate system become: 

Continuity 

" it 7 f pdT ' = f p VV dr ' (21) 

■V 

First Law 




PUp'Eiip'dr' 



du 


Pa dxp' 


dT ' 


(22) 


Momentum 



(23) 


The body force term (J p dT 1 ) g-j-i in Equation (23) is the only term 
which must be added to the equations of motion to make them valid in 
the moving frame . 

In order to preserve the energy conservation property of the 
AFTON 2P equations in the moving coordinate system, a kinematic work 
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term must be subtracted from the surface work term in the total energy 
conservation equation. The total energy conservation equation in the 
moving frame can be derived from the Momentum equation (Equation 19) 
and First Law Equation (20). First the kinetic energy equation in the 
moving frame is derived by multiplying Equation (19) by , i.e.. 

Kinetic Energy Equation 


a 

at- 



+ 3 (p u ' ^ ) 

v~ ) 

+ a Xg ' r u 3 2 / 


dv 


u ■ - ~ 


(24) 


Adding Equations (18) and (24) yields the total energy equation in 
partial differential form 




(E + 


9x a -r^ 


. pV(e + 


U 1 U 1 1 
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dV 
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(25) 


By utilizing the Gauss Theorem, Equation (20), Equation (25) can be 
put into the integral form 


<E + TT' » dT ‘ * f r ' p"b 

-/, V“ r 


* (E +^r~) y dr ' 


dV . 

dT ') dt 


(26) 
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where it can be shown that P-, dr 1 =P D Xo'dT'. Equation (26) is the 

“ pa P 

required integral formulation of total energy conservation in the 

moving coordinate system. It is seen from Equation (26) that the term 
r dV a 

( j ,p U - dT') - dt; must be subtracted from the surface work term 
in order to make the present total energy conservation relation for a 
Eulerian coordinate system valid in the moving frame. 
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AFTON 2P FINITE DIFFERENCE EQUATIONS AT INTERIOR MESH POINTS 
RELATIVE TO A MOVING COORDINATE SYSTEM 

Based on Equations (21), (22), (23), and (26) of Appendix A, the 
AFTON 2P finite difference equations for an Eulerian coordinate system may 
be converted to the moving coordinate system. The field of motion in 
AFTON 2P is actually covered with two closely related finite difference 
meshes -- one for the calculation of thermodynamic variables such as stress, 
i.e., quadrilateral zones, and the other for the calculation of kinematic 
variables like momentum, i.e., momentum zones. 3 Figure 25 illustrates 
the two types of meshes in two space dimensions. The continuity and first 
law equations are applied to calculate properties on a quadrilateral zone 
while the equations of momentum and total energy conservation are used to 
calculate properties on a momentum zone. The finite difference equations 
relative to the moving coordinate system are as follows: 

Definitions : 

R = R(x,y) Cartesian coordinate position relative to the 

moving frame 

U = U(u,v) material velocity relative to a coordinate point 

in the moving frame 

E internal energy 

4 unit vector in the x-direction (i.e., the plane of 

flow, and normal to the free stream flow relative to 
the laboratory frame) 

H total energy 

W rate of work 

V velocity of the moving frame relative to the 
laboratory frame 

p material density 

k unit vector normal to (x,y) plane 

m mass 

V volume 

M momentum 
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M- 

At 


Y 

P 


T ' 

xy 


P 


CT" 


viscosity 
time step 

ratio of specific heats 
fluid pressure 
viscous stresses 



a 

P - CT" 




Ax(t) 

displacement in the x-direction of 
coordinate system relative to the 
(displacement of cylinder relative 
frame) 

the moving 
laboratory frame 
to laboratory 


A 

amplitude of 

cylinder displacement 



U) 

frequency of 

cylinder displacement 



Equations : 




Ax = 

A sin u) t 



(27) 
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II 

A 

1 ui A cos (u t 
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dV 

— c _ 
dt 

A 2 

A A id sin a) t 



(29) 

R = R J 

= R° 



(30) 

U 

U(u,v) 



(31) 

Al 2 = 

(Ri - Rg) x k 



(32) 

A 33 = 

(R 3 - Rg ) x k 



(33) 

Wl3 = 

%(Ui +U 2 ) 



(34) 
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W 33 = %(U 3 + Us) 

(pWA ) 12 = Px 2 W12 • Aig 
(pWA ) 32 = p 3 2 W 3 g • A 32 
V = V 1 = V° 

m 1 = m° + At[-(pWA ) 12 + (pWA ) 32 + (pWA ) 43 - (pWA) 41 ] 

p 1 = m 1 /V 

d ix = x i( y i-l - y i +2 ) + X i-1 ( y i+2 - y i) + X i+2 ( y i “ y i-l) 

where: i=l,2,3,4 and, for example, when i=l, i-l*=4 

d. = -d 

iy ix 

A ix = * [ (y i-l - y i+2 )/d ix + < y i-l ' y i+l> /d i+l X 

+ (y i+2 ' y i-l )/d i+2 x] 

A iy * % [ (x i-l - x i + 2 )/d iy + (x i-l ' x i-l )/d i + l x 

+ (x i+2 ‘ X i-l )/d i+2 y] 

= A lx u l + A 2 x U 2 + A 3x U 3 + A 4x u 4 

= A- u- + A 0 u 0 + A„ U- + A, u, 

oy ly 1 2y 2 3y 3 4y 4 

(3 v 

JT~ = A. v. + A, v, + A, v, + A. v, 
ox lx 1 2x 2 3x 3 4x 4 
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(40) 
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(43) 
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+ 16 [(P WA) 56 + (PWA) 41 + (PWA) 67 + (pWA) l 8 ] % + V* 

+ ^ [(pWA ) 41 + (pWA ) 32 + (pWA) 18 + (PWA) 29 ] (U 1 + U 2 ) 2 

W =-[% (U 1 + U 6 ) - (F b b+Z bc >+% (Di+V^cc+y 

+ % (U, + % (Mi+V^aa + Vl 

. d I 

H 1 = H° - WAt - ^ • nip UAt - (Hg + V) At 


(73) 


(74) 
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APPENDIX C 


DOWNSTREAM AND LATERAL BOUNDARY CONDITIONS 


The downstream and lateral characteristic boundary conditions used 
in this program were based on two priricipal assumptions, namely, that each 
streamline of the flow field in the neighborhood of the boundary at any 
instant of time constitutes a region of transient one-dimensional slab flow, 
and that the sound speed has its free stream value at the downstream boundary. 
For the case of the moving cylinder, this characteristic boundary condition 
was applied relative to a laboratory frame. The flow field variables can then 
be calculated at either the downstream or lateral boundaries as follows: 

First, knowing the velocity field at an instant of time, the position 
of any particle can be updated by one timestep from an "earlier time" to a 
"later time". In particular, a particle position can be found at an earlier 
time such that one timestep later the particle will arrive at a given boundary 
mesh point. For this purpose the discrete velocity field is assumed constant 
during the timestep and is made spatially continuous by interpolation; we 
interpolated linearly. Thus, one can calculate the direction of the streamline 
along which that particle travels which arrives at the given boundary point 
at the later time. Let the boundary mesh point be denoted as point "B", and 
the particle position at the earlier time as point "P"; also, let the super- 
script "o" define a property at the earlier time, let the superscript "1" 
define a property at the later time, and let the absence of a superscript 
indicate a property defined at a time halfway between the earlier and the later 
times. Due to the transverse motion of the cylinder^ the x-component u° of the 
velocity of the moving coordinate system is given by 


o 

u 

c 


U3A CosUJt 


( 75 ) 


where A is the amplitude of the oscillation and cu is the frequency. The 
required streamline direction can be specified in terms of the components 
u, v of the average velocity vector ^relative to the laboratory frame j 
between points B and P,as follows: 


1 + (3v/dy)°(At/2) (u° + u°| -(3u/dy)° v° B (At/2) 


/ g 


( 76 ) 


1 + (3u/3x)°(At/2)v°^ -(3v/3x)°(u° + u°) (At/2) 


/g 


(77) 


g = 


ft 


1 + (du/3x)°(At/2) 1 + (dv/3y) °(At/ 2) 


- (dv/dx)°(du/3y)°(At/2) 2 (78) 


where u is the average velocity component in the x-direction, v is the 
average velocity component in the y-direction, u° is the x-component of 
velocity at point B, v is the y-component of velocity at point B, At is 

B 

the time increment, and (3u/3x) , (3v/3x) , (9u/dy) are the velocity 


derivatives 
cosines of 

evaluated at the zone centroid (see Appendix B) . 
the required streamline are then: 

The direction 


X DIR 

= u/U 

(79) 


y dir 

= v/U 

(80) 

where 

u = 

\A 2 + ^ 2 

(81) 


Assuming that sound signals travel along the streamline as in one- 
dimensional slab flow, it is also possible to calculate from the known 
earlier-time flow field both the particle velocity and the sound speed at 
a point S from which a sound signal would have to depart at the earlier 
time in order to reach the boundary point at the later time. Based on the 
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zone sound speed (which is assumed uniform in each zone) and the average 
velocity vector between points P and B, the location of the point S can be 
determined approximately from the following relations: 


(x s - x B ) = -(U + C) X DIR At 


(y s - y B ) = -(U + C) y dir At 


( 82 ) 

(83) 


where C is the zone sound speed. From the Equations (82), (83) and the 
discrete velocity field at the earlier time (i.e., time "o") the components 
of the velocity vector at point S at the earlier time are 


u° = u° + u° - (U + C) At 


X DiR (3u/dx) ° + Y nii (Su/dy) ° 


[’ 


v° = V° - (U +C)At(x DIR (3v/9x)° + Y t , td Ov/3y)° 


DIR 


(84) 

(85) 


”s° ■ vKr (' 


r O , O 

|U S + u c 


( 86 ) 


where u° is the x-component of the velocity at S, v° is the y-component of 
b b 

the velocity at S, and Ug is the magnitude of the velocity vector at S. 

Based on the assumption that the sound speed has its free stream 
value at the downstream boundary point B, the sound speed at the point S can 
be determined at the earlier time by linear interpolation. For the downstream 
boundary, the relation is as follows: 


C- = C ro + 2 (U + C) Y MR At^-^ 


(87) 


At a lateral boundary Equation (87) is replaced by the equation 


C s = C ffl + 2 (U + C) X MR At 


£ ~ 

- x. 

B A 


( 88 ) 
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where C° is the sound speed at S, C^,, is the free stream sound speed, C is 
the sound speed of the zone containing the point S, and y^ is the y-coordinate 
of the interior mesh point "A" adjacent to the boundary point B. 

According to the method of characteristics as applied to the one- 
dimensional linear isentropic flow of a polytropic gas, the Riemann variable 

R = U + 2C/(y-l) (89) 


is constant along any sound signal trajectory, where Y is the ratio of heat 
capacities for the gas. By hypothesis here, each streamline in the neighbor- 
hood of the downstream boundary is a region of one-dimensional, linear isen- 
tropic flow during a timestep of the numerical calculation. Hence, the Riemann 
variable R has the same value at the point S at time o, as it has at the point B 
at time 1, namely 

R = U° + 2C°/(Y-1) (90) 

From the assumed free stream speed of sound at the boundary point B, and 
the known Riemann invariant, the particle velocity in the streamline direction 
can then be computed at the point B at the later time, as follows: 

Ug = U° + (C° - C OT )/(Y-1) (91) 


Assuming the direction of flow to be the same at times o and 1, the particle 
velocity can be found as a complete vector quantity according to the equations 


u 


1 

B 


U B X DIR 


(92) 



U 1 Y 
B DIR 


(93) 
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To complete this description of downstream boundary flow, mass and 
internal energy are transported across the boundary at the densities charac- 
teristic of the interior of the zone containing the point S. This boundary 
condition has now been applied downstream in many calculations of viscous 
compressible flow around obstacles, where it appears to provide a good approxi- 
mation to flow out of the region of calculation. 
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APPENDIX D 


NOMENCLATURE 


Cylinder radius 


Amplitude of oscillation of the cylinder 
Local speed of sound, c = 


Free stream speed of sound, c = 


YP_ 


C D’ 


S’ 


C L’ 


V 


Instantaneous drag coefficient, C. 


D 


Dra S 




- i fi 

Time-averaged drag coefficient, = — I CL dt 

D TJ Q D 


Instantaneous lift coefficient, C = 


Lift 


L %p m U^ d 

a <» CO 


- 1 XT 

Time-averaged lift coefficient, C = — I C dt 

L 1 L 


C (r.m.s.), Root mean square lift coefficient, C (r.m.s.) 




(r.m.s. ), Root mean square lift coefficient for a stationary cylinder 
o 


V 


p - p 

a 

Local pressure coefficient, C = "! TT s 

P %Poo U a, 


V 


C , 
v 


Specific heat at constant pressure 
Specific heat at constant volume 
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d, 

M, 

M 

CC 

n, 

P, 

P 

CO 

p, 

R, 

s, 

S, 
s 

o 

t, 

T, 
T, 

u. 


U, 


Cylinder diameter 
Mach number, M = U/a 

U /c 

CO < 


a + U 


Strouhal number for stationary cylinder 
Time 

Temperature 

Time integration interval for evaluation of time averages 

Local velocity component in the boundary layer parallel 
to the wall 

Magnitude of free stream velocity vector 
Magnitude of local velocity vector 


Free stream Mach number, M = 


Shedding frequency 


Pressure 


Free stream pressure 


Riemann invariant, P = 


Y - 1 


Reynolds number, R 


pUd 


Arc length 


Strouhal number, S = 


nd 
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Local velocity vector 


U, 

x. Coordinate normal to free stream velocity vector 

y, Coordinate parallel to free stream velocity vector 

Y, Ratio of specific heats, Y = C^/C^ 

6, Boundary layer thickness 

p. Density 

0 , Dimensionless cylinder oscillation frequency, 0 = 

p^, * Free stream density 

6, Angle measured along the circumference of the cylinder from 

the forward stagnation point. 

T, Characteristic time; the number of cylinder diameters that 

a particle traveling with free stream velocity would move 
in a given time. 

CD , Angular frequency of cylinder oscillation 
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10.00 -8.00 -&C0 -4.00 -2.00 0.00 2.00 4.00 E.00 8.00 10.00 



Figure 1. Finite difference mesh (43 x 80); This 

finite difference mesh has an upstream bound- 
ary located at y =-20 cm, although only the 
portion of the mesh downstream of y =-7 cm 
is indicated. 
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Wall shear stress distributions vs. time at points on the cylindrical surface, plus and 
minus 90 degrees from the forward stagnation point of the cylinder; Problem 133.0, 

R = 100.0, Mg, = 0.20; time is measured from the start of cylinder motion, i.e., 1.19 msec 
after the start of motion. 
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Velocity Vector Plot of the Moving Cylinder 
Flow Field at a time when the free stream 
particles have moved 7.08 cylinder 
diameters; Problem 133.0, R = 100, = 0.20; 

Oscillation frequency 2,42 radius/msec. Os- 
cillation amplitude 10% of cylinder radius. 
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Figure 5. Velocity Vector Plot of the Moving 
Cylinder Flow Field at a time when 
the free stream particles have moved 
7.95 cylinder diameters; Problem 133.0, 

R = 100, M w = 0.20; Oscillation frequency 
2.42 radius/msec. Oscillation amplitude 
10% of cylinder radius. 
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Figure 6. Velocity Vector Plot of the Moving 
Cylinder Flow Field at a time when 
the free stream particles have moved 
8.90 cylinder diameters; Problem 133.0, 

R = 100, = 0.20; Oscillation frequency 

2.42 radius/msec, Oscillation amplitude 
10% of cylinder radius. 
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Cylinder Flow Field at a time when the 
free stream particles have moved 10.45 
cylinder diameters; Problem 133.0, 

R = 100, = 0.20; Oscillation frequency 

2.42 radius/msec, Oscillation amplitude 
10% of cylinder radius. 



Figure 8 

Trajectories of the Vortex centers of 
Problem 133.0; R = 100, = .20; cylinder 

oscillation frequency corresponds to a 
Strouhal number of .17; the double amplitude 
is 10% of the cylinder diameter. 
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Figure 9 

Verity Vector Hot of the moving cylinder 
flow field at a tune when the free stream 
parades have moved 7.17 cylinder diameters; 
Problem 133.2, R = 100.0, M„ = .20; oscilla- 
tion frequency 10 radians /msec, oscillation 
amplitude 107. of cylinder radius. 
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Figure 10 

Velocity Vector Plot of the moving cylinder 
flow field at a time when the free stream 
particles have moved 7.65 cylinder diameters;- 
Problem 133.2, R = 100.0, = .20; oscilla- 

tion frequency 10.0 radians/msec, oscillation 
amplitude 10% of cylinder radius. 
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Figure 11 
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Velocity Vector Plot of the moving cylinder 
flow field at a time when the free stream ^ 
particles have moved 8.925 cylinder diameters; 
Problem 133.2, R = 100.0, M- = .20; oscilla- 
tion frequency 10.0 radians/msec, oscillation >■ 
amplitude 10% of cylinder radius. 
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Figure 12 

Velocity Vector Plot of the moving cylinder 
flow field at a time when the free stream 
particles have moved 10.65 cylinder diameters; 
Problem 133.2, R = 100.0, =.20; oscillation 

frequency 10.0 radians /msec , oscillation ampli- 
tude 10% of cylinder radius. 
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Figure 13 

Wall shear stress distributions, tis. time at 
points on the cylindrical surface plus and 
minus 90 degrees from the forward stagnation 
point; Problem 133.2, R = 100.0, M = .20; 
time is measured from the start of cylinder 
motion, i.e., 1.19 msecs after the start of 
motion. 
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Figure 15 

Strouhal number variation with dimensionless 
cylinder oscillation frequency; R = 100.0, 

M =0.2; double oscillation amplitude equal 
to 10% of the cylinder diameter. 
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Figure 17 

Pressure coefficient distribution along the 
cylindrical surface at a time when the free 
stream particles have moved 10.45 cylinder dia- 
meters; Problem 133.0, Instantaneous Lift 
Coefficient, = 0.1665. 
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o AFTON 2P Numerical Data, Oscillation Amplitude 
107 o of Cylinder Diameter 

□ Cincotta's Experimental Data, Oscillation 
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Figure 19 

Lift amplification factor as a function of the 
ratio of the cylinder excitation frequency to 
the shedding frequency for a stationary cylinder; 
o. data points indicate AFTON 2P calculations at a 
Reynolds number of 100.0 and a Mach number of .20. 
□ data points indicate experimental results at a 
Reynolds number of 6xl0 6 and at a Mach number 
of .30. 
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xs .. Figure 20 

Velocity Vector Plot of the moving cylinder 
flow field at a time when the free stream 
"“ " particles have moved 5.08 cylinder diameters; 

Problem 133.1, R = 1000.0, M^ = .20; oscilla- 
. v '.^:“ tion frequency 2.42 radians/msec, oscillation 
amplitude 10% of cylinder radius. 
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Figure 21 

Velocity Vector Plot of the moving cylinder 
flow field at a time when the free stream 
particles have moved 7.625 cylinder dia- 
meters; Problem 133.1, R = 1000.0, M = .20; 
oscillation frequency 2.42 radians /msec, 
oscillation amplitude 10% of cylinder 
radius. 
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Figure 23 

Velocity Vector Plot of the moving cylinder a 
flow field at a time when the free stream > 
particles have moved 9.68 cylinder dia- 
meters; Problem 133.1, R = 1000.0, 

= .20; oscillation frequency 2.42 radians/ 
msec, oscillation amplitude 107. of cylinder ^ 
radius. 
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Figure 25 

Schematic Diagram of an AFTON 2A or 2P Space 
Mesh, Showing a Momentum Zone and the Four 
Quadrilateral Zones Whose Corner Pieces Make 
Up the Momentum Zone. 
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